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Abstract 


The present work deals with the the numerical simulation of turbulent 
flow using a finite volume method and a standard k — e model. A standard k — e 
model was incorporated in the laminar code and was applied on two test cases 
1) Channel flow 2) Backward facing step flow. The code is three dimensional 
and can simulate a flow in a any irregular geometries. The equations are solved 
in a physical domain thus avoiding a need for transformation of governing 
equations. 

The accuracy of determining the turbulent flow field depends upon the 
type of numerical scheme used. The code is dveloped for 3-D applications but 
the present computaions were done using the 2-D modification. The present 
work has demonstrated the ability of the method used to predict the turbulent 
flow field correctly. Standard k — e model is seen to establish good compromise 
between accuracy and computational cost. The results of flow over a backward- 
facing step are seen to be in good aggrement with the experimental as well 
as other advance turbulence model results. In fact it was observed that the 
present results was able to capture the secondary recirculation zone which 
previous computation with the same model failed. The length of primary 
recirculation zone obtained is 6.0 H which is found to be in good aggrement 
with the experimental results of Kim et al 37 for which the length is 7.0 H. 
Also the secondary recirculation length obtained is 0.83 H which also aggres 
well with the length obtained by Lien and Leschziner 40 which is 1.0# using a 
Reynolds stress model. 
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Chapter 1 


Introduction 


During the early days of the computational fluid dynamics(CFD) finite- 
difference methods see e.g Patankar 1 were the most popular. They are algo- 
rithimically simple, efficient, and accurate. However they are best used on 
uniform grids and hence in regular computational domains. However, with the 
advancement of C.F.D and its application to industrial problems there is a need 
for methods for computing flows in complex geometries. To adapt the finite 
difference method for complex geometries, we can map the complex domain 
onto simple domains and solve the problems there. This requires that the gov- 
erning equations be represented using generalised curvilinear co-ordinates this 
makes the Partial differential equation quite complicated and leads to a loss 
of computational efficiency and accuracy. Alternatively one can use schemes 
based on the Finite element or Finite volume methods. 

Finite element methods have been increasingly applied to CFD see e.g 
Hughes 2 . They can be applied to complex flow domains. However these meth- 
ods though very powerful are not simple to implement and their computational 
efficiency vis-a-vis the finite difference method is doubtful 
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An alternate development is that of the finite volume methods see e.g. 
Hirsch 3 which are essentially a generalisation of the finite difference method 
but uses the integral form of the governing equations of flow, rather than their 
differential form. This gives greater flexibilty in handling complex domains as 
the finite volumes need not be regular (rectangular) shaped. This approach, 
too, has some problems the staggered grid with non-coinciding pressure and 
velocity grid points, traditionally used to by finite difference solver to avoid nu- 
merical pressure-velocity decoupling are difficult, and computationally expen- 
sive, to use in complex finite volume codes. Thus the use of the non-staggered 
or collocated grid, with the same locations for pressure and velocity grid-points 
is becoming popular. 

The simplest approach to tackle complex domains is to use an orthogonal 
grid but make modifications on the boundary to accomodate the irregularity. 
The advantage of this is that orthogonal grids are used, but with the disad- 
vantage that the boundaries are only approximately represented. Orthogonal 
grids have often been used for solving the flow equations for a complex ge- 
ometry. Pope 4 applied the conservative forms of the equations with the finite 
volume method for simulating turbulent flow in diffusers. Rapley 5 also used or- 
thogonal grids for calculation of turbulent flows. The calculations in entrance 
region of a chanel with an arbitrary shaped cross section were performed by 
Whitlaw et. al. 6 . 

A general curvilinear co-ordinates are usually used along with non-orthogonal 
body-fitted grids which exactly coincide with the boundaries of the complex 
physical domain. The governing partial differential equations transformed in 
terms of these curvilinear co-ordinates are then solved. The successful use of 
this method was made among others by Galchen and Somerville 7 and Faghri 
et. al. s . Faghri 8 used algebraic co-ordinate transformation. Acharya and 
Patankar 9 used an adaptive grid procedure for parabolic flow problems for 
irregular geometries. All these methods employed staggered grids. The sue- 
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cessful application of finite volume method using a collocated i.e non-staggered, 
grid was made by Hsu 10 ,Rhie u ,Peric 12 and others. Hsu 10 and Peric 12 solved 
the equations on a physical domain thus avoiding the need for transforma- 
tion, while Rhie mapped the non-orthogonal grid onto a rectangular mesh on 
which the transformed equations were solved. Mukhopadhay et. al. 13 used 
the staggared grid for the solution of equations on the physical domain without 
transformtion. Verma and Eswaran 14 have developed a Finite volume approach 
directly applied on the physical domain. Satya prakash et. al. 15 developed 
a transient Navier-Stokes algorithm using finite volume method employing a 
collocated grid for complex geometries, in which the equations are solved in a 
physial domain. 


1.1 k-e modelling of turbulent flows 


The numerical solution of the transient Navier-Stokes equations for turbulent 
flows without turbulence modelling , but with sufficient accuracy and reso- 
lution to capture turbulence behaviour, are called Direct Numerical Simula- 
tions(DNS) of turbulence. Direct Numerical Simulation, although an unsur- 
passed method for understanding the physics of turbulent flow can be currently 
used only for Reynolds number less than 10 3 due to constraints of computer 
memory and speed. Practical needs of industry demands numerical solutions 
for high Reynolds number (10 6 ) flow problems. 

An alternative to Direct Numerical Simulation is to use the Reynolds 
Averaged Navier-Stokes equations. In this approach, the statistical average 
of the Navier Stokes equations are used as governing equations to describe 
the mean quantities. However, the averaged equations no longer constitute 
a closed system of equations. The system can be closed only with the aid 
of semi-emperical models. Averaging gives rise to extra stress terms called 
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Reynold stresses, which have to be so modelled. 

Usually in turbulence modelling, the Reynolds stress terms are approxi- 
mated using a Boussinesq eddy viscosity concept which assumes that, in anal- 
ogy to the viscous stresses in laminar flows, the turbulent stresses are propor- 
tional to the mean velocity gradient. The whole purpose of the modelling then 
becomes the prediction of the eddy viscosity which substitute the kinematic 
viscosity of the Navier-Stokes equation. The simplest model for eddy viscos- 
ity is the Constant Eddy Viscosity model which however has few successful 
applications. Another method is to specify an algebraic relationship for eddy 
viscosity as in e.g., the Prandtl mixing length model. Such models are called 
Zero equation models because no additional partial differential equations are 
needed to be solved to implement the model. 

One equation models specify the eddy viscosity by an emperical rela- 
tionship with turbulent kinetic energy and a charecteristic length-scale. The 
turbulent kinetic energy is computed using a partial differential equation and 
the length-scale is specified algebraicaly. 

The most succesful popular model , however, k-e model which is a two 
equation model that uses transport equations for k and e and with these two 
equations obtaine the eddy viscosity. The k-e model contains five constants 
which are not truly universal and thus varies according to the flow situation. 
Rodi 16 lists several cases where the k-e model has been successfully applied 
using given model constants values. 

Some of the difficulties of k-e model are. 

1. The k-e model uses a Boussinesq eddy viscosity concept which assumes 
that turbulent viscosity is isotropic. Therefore it cannot reproduce sec- 
ondary flows which arise due to unequal normal Reynolds stresses. 
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2. Boussinesq eddy viscocity concept cannot accurately simulate the effect, 
on turbulence, of the changes in the mean flow direction across the shear 
layers. 

3. The turbulence production in highly strained flows is over-predicted by 
the k-e model. 

4. The k-e model requires modification when applied to flows with highly 
curved stream lines. 

5. The k-e model does not strictly satisfy the realizability constraint, (r^ > 

0 ). 

A variety of modified version’s of the k-e model have been used to im- 
prove its performance, among which are the RNG k-e model of Yakhot et. 
al. 17 and the Non-linear k -e model of Speziale 18 . Some Low Reynolds number 
form of k -e model have been developed which use damping functions to damp 
the turbulence near a solid wall surface due to molecular viscosity. Many such 
model’s appeared with different constants k, damping function’s the compari- 
son of which is given in Patel 19 . 


1.2 Development of Finite Volume Code for 
complex 3-D geometries 


A Finite Volume Code to solve laminar Navier-Stokes equations has been de- 
veloped by Satya Prakash et. al . 15 . The code has been tested on different 
regular geometries like l)Lid driven 2-D and 3-D cavity flow. 2)Channel flow 
3) Sudden expansion in a channel. The results in these studies are in good 
agrrement with published results. 

T his code uses an explicit transient algorithm by which both steady and 
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unsteady flows can be solved. The Navier-Stokes equations are solved on a 
non-staggered, and possibly non-orthogonal, grid . A special (momentum) in- 
terpolation technique is used to avoid the decoupling of pressure from the 
velocity field in the numerical solutions. An implicit time stepping version 
has also been developed. Efforts are in progress for testing this code on a 
cylindrical and irregular geometries. 

In this work presented here we incorporate the standard linear k-e model 
of Launder et. al. 2Q with wall function approach of Launder and Spalding 21 
for simulation of turbulent flows in the above solver. The extended code is 
then tested on channel flow and flow over a backward facing step, results. 



Chapter 2 


Finite - Volume Formulation 


In this chapter we describe the finite-volume formulation used in the 
DRDL code (Satya Prakash et a/ 15 ). The formulation included only the Navier- 
Stokes and continuity equation. The k — e model has been incorporated as part 
of this thesis work. 


2.1 Governing Equations 


The three-dimensional Reynolds-averaged Navier-Stokes equations for incom- 
pressible turbulent flow for an arbitrary spatial region of control volume V 
bounded by a closed surface S can be expressed in the following general 
convection-diffusion-source integral form: 

I u • dS = 0 (2.1) 

Js 

| X ^ + Jj# - 7 V *] ■ ds = U ^ (2 - 2) 

where p is the fluid density, u is the mean fluid velocity with components 
it, v, w, (j> stands for any vector component or scalar quantity, S# is the volu- 
metric source term. The various source terms are. 



2.2 Description of the Finite- Volume 


8 


s w 

Sk 

S< 


1 dp d , du . d . dv . d . dw. 

~pdi + + + ai^dJ 

1 dp d . du. d , dv. d . dw. 

~~pd y + dx^dy^ + dy^dy) + 'dz yt ~dy ) 
1 dp d , du. d . dv. d . dw. 

- pTz + + di^dJ + 


CurG - cjr; 

k k 


( 2 . 3 ) 

( 2 . 4 ) 

( 2 . 5 ) 

( 2 . 6 ) 
( 2 . 7 ) 


In this formulation we work with cartesian components of velocity. So <j> can 
be the three cartesian components of velocity u, v, w, as well as any scalar e.g., 
k, or e which needs to be determined. The last two source terms correspond to 
the equations for kinetic energy k ,and the turbulent energy dissipation e,from 
the k-e model, C\ t and Ci t are model constants and G is the generation term 
which shall be explained later. 


2.2 Description of the Finite- Volume 


Figure 2.1 shows the kind of control volume used, with eight vertices but 
otherwise irregular. 

The conservation equations are discretized by employing the general 
finite-volume approach. The solution domain is subdivided into a number 
of contiguous (finite) control volumes. The control volumes are defined by the 
coordinates of their vertices which are assumed to be connected by straight 
lines. The coordinates of the control volume vertices are calculated by some 
grid generation procedure. A Collocated Grid arrangement is employed, i.e., 
all the dependent variables ( u , v,u>,p, k, e) are defined at the same location, the 
centroid of the control volume in Figure 1. The six neighbouring control vol- 
ume centers are indicated by E, W, N, S, T, and B (for the east, west, north, 
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south, top and bottom neighbours). The face center points e, w, n, s, t, and b 
are located at the corresponding face centroids of control volume P. Note the 
edge-centers te, be, ne, se, etc. which are also needed in the computations. 

Before the time-stepping can proceed all geometrical parameters need 
to be computed as part of the initialization. These parameters include the 
surface vectors for each of the six faces e, w, n, s, t, b of each finite volume and 
also its volume. These computations can be done by the methods suugested 
by Kordulla and Vinokur 22 and the details of the application can be found in 
reference 15. 


2.3 Discretization Procedure 

All the conservation equations have the same general form, represented by 

^ J v p<t>dV + J s [pM<j> - r '+V<f>\ • dS = J v SfidV (2.8) 

The main steps of the discretization procedure to calculate convection and 
diffusion fluxes and source terms are outlined below. The rate of change and 
source terms are integrated over the cell volume, whereas the convection and 
diffusion terms form the sum of fluxes through the CV faces. 

2. 3.0.1 Continuity Equation 

Equation (2.1) is discretized in the following way. 

J pu • dS =» 53 p( u • S )j (2.9) 

§ j=e,w,n,a,t,b 

where Sj is the surface vector representing the area of the j tfl cell face and u, 
is the velocity defined at the face center j. 
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In discretized form the continuity equation follows 

Y / F J =F e + F w + F n + F 3 + F t + F b = 0 (2.10) 

i 

where the F } is the outward mass flux through face j, defined by 

Fj = pUj • S j 


2.3.0.2 Discretization of General Convection-Diffusion Equation 


(a) Rate of change 

The value of the dependent variable <f> at the centroid of the control volume 
(the geometric center) represents an average over the CV as a whole. Thus, 
using forward difference time-stepping we get: 


d_ 

dt 


J v P<fidV^pV 


At 


( 2 . 11 ) 


where V is the volume of the cell. 


(b) Convection fluxes 

The surface integral of the convection flux of variable (f> is approximated in the 
following way. 

f pu(/)-dS => 53p^(u-S)j (2.12) 

JS j 

= 

j 

where <f>j is the value of 0 at the center of face j. Thus 

/ pu<£.dS => F e <f> e + F w (j) w + F n <f> n + F s <p s + Ft<f>t + F^b (2.13) 

Js 

where <f> e is the (interpolated) value of the variable <j> at the east face center, 
etc. This is evaluated by using a weighted upwind interpolation between the 
neighbouring nodal values, say <frp and 4>e ■ This weighted upwind methods 
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takes weighted average of the centre-difference (linear interpolation) between 
4 ip and 4> E and the first order upwind approximation. 


,uds _ r <t>P if F e > 0 
™ e *■ 4>e otherwise 


(2.14) 


The centre-difference approximation is 

,cds _ Ve4>p Vp<t>E 
9 V E + V P V E + V P 


(2.15) 


where V E and Vp are the volumes of the cells E and P respaectively. Here cell 
volumes are used for interpolation, instead of cell dimensions, so as to handle 
irregular finite- volumes if necessary. In the weighted upwind discretization the 
convection flux is split into the first-order upwind differencing approximation, 
and another part, which equals the difference between the central difference 
scheme (CDS) and UDS approximations: 


F e <t>c = (F e <j> e ) UDS + 7 l(F e <f> e ) CDS - (FM UDS ], etc (2.16) 


Multiplication of the second part by a factor 7(0 < 7 < 1) allows the 
introduction of numerical diffusion. (7 = 0 means pure UDS, 7 = 1 means pure 
CDS) .The lower value of the 7 enhance the numerical stability of the solution 
algorithm, but at the cost of accuracy. In practice one should use the highiest 
value of 7 < 1, which allows numerical stability. This optimum value of the 7 
depends upon the grid peclet number of the computations. If a fully implicit 
method is used for time-stepping, the upwind parts of the equation(2.16) are 
‘implicit’ in that they are incorporated in the coefficients of the unknown 
velocity during the pressure-velocity iterations. The CDS terms on the other 
hand are evaluated using previous iteration values and used as source term 
on the right hand side of the same equation. This is the so-called deferred 
correction approach of Khosla and Rubin 23 which has been incorporated in 
the implicit code(reference 15). 

(c) Diffusion fluxes 
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The diffusion flux of variable 0 through the cell faces is evaluated as given in 
reference 15. 

jfr*v«s-ds* £ (r*v,i • s), = £ -Ff (2.17) 

j~e,w,n,s,t,b j 

where any face S, can be expressed as 

S_, = ain 1 + a 2 n 2 + a^n 3 (2.18) 

where n 1 , n 2 and n 3 are any three linearly independent (not necessarily orthog- 
onal) unit vectors. The details of the calculations of ai, 0:3 which depends 
on the local grid perticulars can be found in Satya Prakash et al. 15 . 


The diffusion flux is made of two distinct parts: normal diffusion and 
cross-derivative diffusion. The second part arises from the nonorthogonality of 
the grid. The normal derivative diffusion flux approximation of 0 through any 
cell face involves the values of 0, at cell centres whereas the cross-derivative 
diffusion flux approximation takes into account the edge center values of 0. 
(The normal derivative diffusion flux is treated implicitly in implicit time- 
stepping while the cross-derivative diffusion flux is treated explicitly to avoid 
the possibility of producing negative coefficients in an implicit treatment). The 
diffusion flux at any face of the cell, say east, is given by. 


F e — — r^(ai 


0 E 


Ax 1 


0P , 0se 
4- 0:2 


Ax 2 


0ne . 0te 06c \ 
+ q 3 — r ~) 


Ax 3 


(2.19) 


The edge center values appearing in cross-derivative diffusion flux are 
calculated by using the linear interpolation (Satya Prakash et. al. 15 ). 

(d) Sources 

The source term is integrated over the cell volume. By assuming that the value 
at the CV center represents the mean value over the whole control volume, we 
can write 

J^StdV => (S^pV 


( 2 . 20 ) 
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Apart from the real source S$, explicitly treated parts of the convection and 
diffusion fluxes are added to S^. 

The pressure term in momentum equation is treated explicitly in the 
predicted step(see below). Its discretization for the i th momentum equation is 

~ JgP 71 ' S =^-Y^PjSij 

where p 3 is the pressure at the j th face center and S t j is the i th component of 
the surface vector for face j. 

2.4 Semi- Explicit Time Stepping 

We use an explicit time-stepping algorithm for the scalars A, e (and any others 
we may want to compute later temperature etc.). That is every term for the 
convection, diffusion and source are evaluated using current known values at 
next time step V to compute the values at the next time step, explicitly, as this 
unknown value appears only once in the time-derivative. This simple approach 
cannot be adopted for the momentum equations because the pressure and the 
velocity components are simultaneously unknown. The pressure has to be 
such that the continuity equation is satisfied, but there is no seperate equation 
for pressure. Thus iteration is inevitable and the approach adopted is semi- 
explicit-i.e. the convective, diffusive and source terms are treated explicitly i.e. 
evaluated at ‘n’ time level, but the pressure is treated as unknown. 

For the present situation we will adopt a semi-implicit scheme in which 
the three momentum equations 

pv ^ - ~ + UF‘ + < 2 ' 21 ) 

J j 

where u is the i th component of velocity, have to be solved along with 

£ F? +1 = 0 


( 2 . 22 ) 
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for each finite volume cell. We adopt a two step process. First a predicted 
velocity u* is found which satisfies the equation 

pv ^zr + ? (J?C + p- 23 ) 

Subtracting equation (2.23) from (2.21) we get 

< 2 - 24 ) 

where the corrections 

u'p = Up +1 — Up Pj~ P" +1 — P'j (2.25) 

and the corresponding flux corrections Fj have to satisfy 

= ( 2 - 26 > 

3 j 

In the predictor step, the predicted velocity is computed using(2.23),and in 
the corrector step(2.24) and (2.26) are treated till convergence is achieved and 
and the correction of pressure and velocity are added to the predicted values 
to obtain the final values for that time-step. 

u n+1 = u* + u' (2.27) 

p n+ 1 = p n + p' (2.28) 

Due to the non-staggered variable arrangement, if the variables (velocities 
and pressure) at the cell faces are calculated by linear interpolation between 
the adjacent cell centered quantities then the pressure-velocity iterations do 
not converge and lead to a checker board pressure field. Therefore it is impor- 
tant to use momentum interpolation(HiTshet. al . 24 in which the velocity at the 
cell faces are computed by allowing linear interpolation of the convective and 
diffusive terms but not of the pressure term. 
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2.5 Solution Algorithm 


The velocity and pressure fields are calculated with the following Gauss-Seidel 
type algorithm: 

1. Use equation (2.23) to compute the cell-center up and Vp. Make an 
initial guess pp. Use interpolation to obtain the face-center quantities v, 
and Pj. 

2. Compute the mass flux through each cell face j using 


Ff = prvj • Sj — A tVpj ■ S j (2.29) 


3. Use equation u'- = — ^Vp' to compute the flux correction at the face j, 

F; = —A tVp'j ■ Sj (2.30) 

This is computed using the formulation for diffusion fluxes by replacing 
<t> by p'. 

4. Compute residual for each cell, 


(2-3i) 

3 3 

5. Calculate the cell-center pressure correction from the relation 

Pp = Pp + oj— (2.32) 

dp 

where u is relaxation factor and ap stands for diagonal coefficients that 
is calculated from 


(2.33) 

where oq, et 2 , etc. are the same as in p gradient calculations above. 


ap = —At 


ai 


c*i 


LAa : 1 


+ 


a 2 


a 2 


Aa : 1 e Aa; 21 " Aar 2 " Aa ; 3 


» + 


a 3 


<*3 

Aa : 3 
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6. If > e goto step 3 

7. Store the updated mass flux through cell-faces from 

F* = F*~ A tVp] ■ S 3 (2.34) 

8. Store the updated pressure at cell-center pp — pp + p'p 

9. Store the cell-center corrected velocity 

Up = (2-35) 

Up = up 4- Up (2.36) 

10. Calculate the turbulent viscosity by solving the k and e equation. 

11. Goto step 1 and repeat the process until steady state is reached. 



Chapter 3 


k — e Modelling 


The continuity equation and the Navier-Stokes equation, with the asso- 
ciated boundary conditions, completely describe the behaviour of a Newtonian 
viscous fluid flow. For an incompressible fluid, the governing equations are 


continuity equation 



(3.1) 


and the Navier-Stokes equation 

du t dui _ 1 dp dPui 

dt + Uj dxj pdx x dxj dxj ’ 


(3.2) 


where u t is the instantaneous velocity component in the x t direction, p is 
the pressure, p is the density, v is the kinematic viscosity, and the Einstein 
summation convention applies to the repeated indices (i and j take a value of 
1 to 3). Although they are used mainly to solve laminar flows, equation (3.1) 
and (3.2) are valid for any Reynolds number, i.e., even for turbulent flows, but 
in this case the solution fluctuates so rapidly, and at such small scales, that it 
is extremely difficult to obtain an accurate numerical solution. 
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However such solutions (called Direct Numerical Simulations) of the 
above two equations can be attempted for turbulent flows, but the difficulty 
is that the capacity of the present computers does not allow the DNS to be 
used for high Re flows. Also engineers are usually not as interested in the 
fluctuating velocity of turbulence as in the mean flow behaviour. Therefore 
a statistical approach is usually taken. This approach usually involves the 
Reynolds Averaged Navier-Stokes equations. 

3.1 Reynolds Averaged Navier-Stokes equation 

To obtain these equations, a time average is taken of the Navier-Stokes and 
continuity equation. Velocities and pressure fluctuation are decomposed into 
mean and fluctuating quantities 

Ui = u + u'i , p = p + p (3.3) 

in which the overbar denotes the mean quantity and the prime denotes the 

fluctuating part. The mean quantites, u \ and p are defined as 

1 r t <>+ T 1 f t o +T 

u, = - / u x dt, Pi = 7f Pdt (3.4) 

1 Jto i Jto 

where t 0 is any arbitrary time and T is a time scale larger than the largest 
time scale of turbulence. During averaging the following rules are valid 

u'i 

uTuJ 

U x Uj 

The averaging then leads to 


= 0 , p '=0 

= U x Uj + u'jUj 


= UjU x = 0 


du x 


= 0 


(3.5) 


( 3 . 6 ) 
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The equations for the average quantities are very s imil ar to the equations 
(3.1) and (3.2) for the instantaneous velocity and pressure fields. The extra 
terms (ujUj) in the above equation, in comparison to equation (3.2) arise due to 
averaging of the nonlinear convective terms (u- Vu) in the momentum equation 
and are known as Reynolds stresses. These terms form a symmetric tensor, 
and are often called stresses because they correspond to the retarding effect 
due to turbulent transport of momentum. 


3.2 Boussinesq Approximation 


The Reynolds-stress tensor (for incompressible flows) is often modelled by the 
widely-used Boussinesq approximation in most present day turbulence models. 
This approximation involves the so-called Eddy- Viscosity concept, and is based 
on an analogy with the viscous stresses in laminar flows, the turbulent stresses 
are also assumed to be proportional to the mean velocity gradients in the 
following way 


U t Uj = -v t 


' dui du. 
. dx a "** dx 


~] + 3 kS,i ’ 


(3.8) 


in which v t is the turbulent or eddy-viscosity. This eddy viscosity v t is not a 
fliud property, like the molecular viscosity v, but is a function of the turbulent 
velocity field. In analogy to the kinetic theory of gases, the turbulent viscosity 
may be expressed as 


v t a v 0 l 0 or v t oc ll/r 0 


(3.9) 


in which no is the relevant velocity scale, la is the length scale and To is the 
time scale of turbulence. However, it is well known that the turbulent flows 
contains a wide range of length and time scales and so there is an inherent 
limitation to the applicability of the eddy viscosity concept. There are also 
other reasons why this concept is criticised, (see, e.g, Biswas and Eswaran 25 ). 


Inspite of these limitations, this concept is very widely used in turbulence 




3.3 Types Of Turbulence Models 


21 


models today. Its great advantage lies in that the problem of modelling the 
six independent quantities, is replaced by the problem of modelling the 
single quantity u t , the eddy-viscosity. Thus specification of the eddy-viscosity 
constitutes the core of the majority of the turbulence models for engineering 
applications. These models are discussed in the next section. 


3.3 Types Of Turbulence Models 


The following models all assume the Boussinesq approximation(3.8). They 
differ only in the method used to estimate the eddy- viscosity, v t . 

3.3.1 Zero-Equation Models 

These are the simplest among the turbulence models. They do not employ the 
use of any extra transport equation for, say, the turbulent kinetic energy or 
some other quantity, for the purpose of modelling. 

Two of the commonly used Zero-equation models are 

(a) Constant Eddy Viscosity Models: 

In these models, v t is specified as a constant and its value is chosen depending 
on some experimental information and or by trial and error procedure . These 
models are obsolete. 

(b) Mixing Length Models: 

This model was propsed by Prandtl. He postulated that 

1. the turbulent length scale, l 0 is equal to the mixing length, l m and 

2. the velocity scale, vq is equal to the mean velocity gradient times the mixing 
length. 
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when these assumptions are substituted into equation (3.9) we get 




du 

dy 


(3.10) 


where the constants of proportionality has been assumed to be unity. This 
model is closed once the mixing length, l m is specified by some simple emperical 
formulation. 


3.3.2 One-Equation Models 


In these models, the velocity scale, Vo of the turbulence is taken as \/k, where 
k is the turbulent kinetic energy k = The eddy viscosity then can be 

written as. 

u t = Vkl 0 (3.11) 


In one-equation models, l Q is specified emperically while the turbulence kinetic 
energy is determined by solving the model equation: 


dk _ dk 
dt + Uj dxj 


v t 


du t du, 

dui d 

+ 

u t dk 

dxj dx t 

dxj dxj 

cj dx j 


p/2 

-C D ~— + uV 2 k (3.12) 
‘o 


in which cr fc and Cu are emperical constants. For simple flows, the suggested 
value of crfc and Co are 1.0 and 0.08 respectively. 


3.3.3 Two-Equation Models 

Two-equation models are among the most popular for scientific and engineering 
calculations. In these models, two separate transport equations are solved to 
determine the length and velocity scales for the eddy viscosity. 

The transport equation for the velocity scale in the two-equation models 
is usually the k equation [equation(3.12)j as given for the one equation models. 
However, the other transport equation need not necessarily have the length 




3.4 k — e Model 


23 


scale itself as the dependent variable - any combination of k and Z 0 may be 
used as the dependent variable in the second transport equation. Kolmogorov 
proposed a transport equation for the frequency Vk/lo, Rotta proposed an 
equation kl 0 , while Spalding suggested a transport equation for the vorticity , 
fc/Zg. The k — e uses a transport equation for the turbulent dissipation rate, e. 


The transport equation for any of the above quantities can be expressed 
in the general form: 


dz _ dz 
dt + Uj dxj 



/ Vklo dz \ 
V dxj) 


+ C zl -G 

k 


C z 2Z—j— + S 


(3.13) 


in which, z = k m lQ, G is the production term of z, S' = is a secondary source 
term which depends on choice of z, and a z , a z y and cr z2 are emperical constants. 
Past experience indicates that the k — e model works better than the other 
two-equation models based on equation (3.13). The description of the k — e 
model is presented in the next section. 


3.4 k — e Model 

The length and time scales in this model are obtained from the turbulent 
kinetic energy, k, and the turbulent dissipation rate, e, using the dimensional 
arguments 

Zo oc — ^ ■ and To a — ; (3-14) 

e e 

substitution in equation (3.9) gives 

u t = C^- (3.15) 

e 

in which is an emperical constant. The turbulent kinetic energy and the 
dissipation rate are obtained by solving their respective modeled transport 
equations. The transport equations for k in the k — e model require mod- 
eling of only the turbulent transport term whch is done using the gradient 
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approximation. This is given as 


convection 


dk _ dk 
dt + Uj dxj 


production turbulent diffusion 

1 A - " -■ S. > . 


du t 

•sTl 

1 

H 

din 

* 

d_ 

f v t dk\ 

dxj 

dXi 

dxj 

dxj 

dXj J 


dissipation of k molecu lar di ffusion 

'T' + vV 1 2 k (3.16) 


The transport equation for e is 


convection 
de _ de 

dt + U] dxj 


molecular diffusion 


turbulent diffusion of e 

iV 

dx 3 \a e dxj ) 


+ 


production dissiaption 



(3.17) 


Speziale showed that the modelling of the producition term in the e 
equation is strictly valid only when 


1. there is a clear separation in the time scale of turbulent fluctuation and 
the mean velocity field and 

2. the level of anisotropy is not very much 

We use the constant values given by Launder and Spalding 21 C fl = 1.44, 
C (2 — 1-92, Ok = 1.0, o e = 1.30 and = 0.09. The standard k — e model 
comprises of equation (3.15), (3.16) and (3.17). These equation along with 
equations (3.6) and (3.7) form a closed system for determining u t , p, k , c in 
an incompressible turbulent flow. 
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3.5 Non-Dimensionalisation of the Governing 
Equations 


The variables are non-dimensionalised in the following m ann er 


U 

Y 


u 

tv 

y_ 

H' 

k_ 

ur 


V = A w = ± x = - 

Uo’ u 0 ' H 


7- z P- p 


Po 




t = 


H/U 0 


6 n — 


e _v t 

VS/H' ” v 


(3.18) 


where Uq, H,po are some charectristic scale of velocity, length, pressure, respec- 
tively and v is the kinematic viscosity. The equations are then transformed 


to 


• Continuity equation 


dU dV dW 
dX + dY + dZ 


(3.19) 


• Momentum equations 

The compact form of the moementum equation in non-dimensionalised 
manner is 


dUj - dUj 
dt 1 dxXj 



where Reynolds number Re = UqH/ v 


• k equation 


Dk n 

~Dt 


i d 
RedX 
1 d 


+ 


RedZ 


( 1+ ^f) 

( 1+f ^f) 


dkn 

dX 

dkn 

dZ 


+ 


i d 


RedY 
+ G n ~ 


1 + 



CEN 



I# 4 
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where the nondimensional production term is 



• e equation 


De^ 

Dt 


1 d 
RedX 



Re dY 





and the non-dimensionalised eddy viscosity i/ t>n is 

v t>n = CpRe— (3.24) 

equation (3.19 - 3.23) constitue the full set of equations which need to be solved 
with the relation (3.24) to obtain the mean flow filed. 


3.6 Boundary and Initial Condition for Tur- 
bulent Flows 

The partial differential equation (3.19 - 3.23) need appropriate initial and 
boundary conditions to yield unique and physically meaningful solutions. 
Boundary Condition 

The inlet condition for the mean velocity can be specified using any physically 
reasonable profile. For velocity we can prescribe a uniform profile U = U 0 ad 
the inlet value of the velocity, or a fully developed channel velocity profile can 
be specified. One may perhaps specify the 1/7 power law velocity profile at 
the inlet. The specified velocity profile at the inlet depends on the type of flow 
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beign simulated. Typically the transverse velocity components i.e., V and W 
are prescribed as zero at the inlet. As for pressure, the gradient of Presurre 
may be set to zero 

dP 

ax=° < 3 - 25 > 

The turbulent kinetic energy k n and its dissipation rate are calculated from 
the value of turbulent intensity specified at the inlet as follows: 

*m = l-5 (U tn I) 2 , (3.26) 


where I is the turbulent intensity, and 


€m00 


XY 

A 


for Y < (A/x) 
for Y > (A/x) 


(3.27) 

(3.28) 


where x = 0-42 which is known as von Karman constant, A is a constant 
usually taken to be equal to 0.09. and Y is the normal distance from the wall. 


Often the turbulence quantities of the domain are uncertain and one is 
forced to estimate these quantities as best as possible. In general the inlet 
turbulence intensity and characteristic length are a function of flow details 
upstream. A very quiescent flow upstream might yield an inlet turbulence 
intensity of 2 — 3% (i.e. 1=0.02-0.03). If the upstream flow involves flow over 
a rough edges or turning, the inlet value intensity may be as high as 10 — 14% 
. A range of k xn and e,„ are used in the literature but for most modelling 
applications the prediction are found to be insensitive to these values (Luo 26 ) 
and this insensitivity was also observed in the channel flow simulations done 
in this study. 

At the outlet boundary the gradients of the quantities along the direction 

of flow can be set to zero. 

dU_ _ dV _ dW = dk n = de n 
dX ~ dX ~ dX ~ dX ~ dX 


(3.29) 
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The pressure at outlet boundary was set to zero. Details regarding the bound- 
ary condition at the wall for velocity, k and e are discussed in the next section 
which describes the near walltreatment for these computations. 

Initial Condition 

In this sudy, we are interested in only the steady state solutions and times 
stepping is only the false transient means of obtaining the final state. Thus 
the initial conditions are prescribed quite arbitrarily as the steady solution 
doesn’t depend upon them. The initial condition for velocities can be set to 
zero over the entire domain, pressure is also set to zero initialy. For the turbu- 
lent kinetic energy and dissipation rate, the inlet values were taken as initial 
condition. 


3.7 Near Wall Treatment 


The velocity variation near a wall in the a turbulent flow has typical features. 
The velocity distribution shows essential differences in the three different re- 
gions. 

1. the viscous sublayer ( Y + < 11) 

2. the Log-law region (11 < Y + < 70) 

3. the fully turbulent region (F + > 70) 

where Y + is the nondimensional distance to the wall and is expressed as Y + = 

UtY 

V 

In the viscous sublayer the flow is dominated by the viscous effects and 
the turbulent effects are negligible i.e. the contribution from the turbulent 
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friction can be completely neglected compared to laminar friction. In the Log- 
law region both contributions are comparable, whereas in fully turbulent region 
the laminar contribution is negligible compared with turbulent friction. 

For high Reynolds number flows the value of Y + may range above 1000. 
The thickness of the viscous sublayer is very small compared to the width of 
the whole domain. However significant changes in velocity as well as tem- 
perature occur in the viscous sublayer it is very important that this region 
be taken into account while performing the computations. This would mean 
that we must take the computational grid small enough to represent the sub- 
layer with accuracy but require much computer memory and time (as done in 
Low-Reynolds number k — e model). Another alternative is to allow a coarser 
grid points in the sublayer, but to simulate the effects of the sublayer and the 
log-law region using the laws of the wall. This is the approach taken by he 
wall function treatment. 

In the standard k — e model the k and e equation are applicable in the 
fully turbulent region ( Y + > 70). So while using a k — e model usual that the 
first grid point be placed at a distance Y + > 30, and the variables at the first 
grid point are specified by the wall function treatment discussed below. 


(a) Momentum Equation 

In the log-law region the universal velocity distribution may be written as 


or 


i In y + + 5.5 (3.30) 

U r X 


^ = -ln (EY + ) (3.31) 

U T X V ' 


where U T = yj**- is a friction velocity, x the von Karman constant is 0.42 and 
E = 9.0 for smooth walls. Using this velocity profile we can compute the wall 
shear stress from the velocity y at first grid point outside the viscous sublayer. 
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This law is not valid in the region of stagnation and separation but then the 
wall stress has little influence on overall flow and thus the little effect on the 
prediction. 


(b) Kinetic Energy 

The Reynolds stresses in the region 30 < Y + < 100 are nearly constant. In 
the region 30 < Y + < 100, the convection and diffusion of u' t Uj are negligible 
so a local equilibrium prevails. This implies that the production and rate of 
dissipation of turbulent kinetic energy near the wall are equal i.e. 

—dU , . 

~ pU V dY =p€ ( 3 - 32 ) 

The Reynolds stresses in this region are nearly constant at the value t w . 


pu'v = T u 


(3.33) 


from this we can write 


— r w dU 
-UV = — = v t — 
p oY 


Inserting Vt — C^ in the above expression we get 

— r k 2 dU 

* = C “TW 

Substituing the value of obtained from (3.32) we get 


(-uV) = Cji 


fc 2 1 
^ (—u'v lS ) 


or 


k = 


J5_ 

cl /2 


in the nondimensional form equation (3.37) is 


k n 


Ur, n 

cl ! 2 


(3.34) 


(3.35) 


(3.36) 


(3.37) 


(3.38) 


Equation (3.38) is used for calculating the kinetic energy at the near wall point. 
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(c) Dissipation Rate 

In the log-law region, the local equilibrium of production and dissipation of 
turbulent kinetic (3.32) and the constancy of Reynolds stresses yield 


t w dU 


6 p dY 


= U'?* 

T dY 

(3.39) 

From log-law (3.31) we obtain 


dU U T 
dY ~ X Y 

(3.40) 

Finally e becomes 

U* C 3 J 4 k 3 ' 2 
e ” ~ xY r ~ xYp 

(3.41) 

in the nondimensional form equation (3.41) is written as 


€p ’ n XY P X Y P 

(3.42) 

where Y p refers to the normal distnace from the wall . 
represents the value of the e at the near wall point. 

The equation (3.42) 


3.8 Implementation of k — e Model 


In this section we discuss the way in which the k — t model has been imple- 
mented, mainly dealing with the the wall- function treatment and the algorithm 
for obtaining eddy viscosity. 


The Reynolds averaged equation is solved over the entire domain with 
a modification when applied to the near wall point. The no-slip condition is 
imposed on the wall and the velocity at near-wall point is obtained from the 
velocity gradient computed fromm the wall shear stress: 


du uxcyw* 
** dy HEyt) 


(3.43) 
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where p refers to the first grid point near the wall. For the application of 
equation (3.43) it is very important that point p be at a distance of y + > 11.63. 


For turbulent kinetic energy and dissipation rate equation (3.16) and 
(3.17) are solved in a domain far from the viscous sublayer i.e. above y + > 30. 
For the near wall point the values of k and e are obtained through the use of 
equation (3.38) and (3.42) respectively. For a three dimensional turbulent flow 
u T is expressed as 


U r = 


xU 


ln(E 9 +) 

where U is the resultant velocity parallel to the wall 


(3.44) 


Solution Algorithm 


1. Guess the initial velocity, k , e and p. Initial values of k and e has to be 
taken as given in section (1.6). The initial turbulent viscosity has to be 
computed from the initial values of k and e. 

2. Solve the momentum equations for the velocity at (n+1) time level using 
the value of u t at (n) time level. 

3. calculate the convection and diffusion fluxes for k and e as described in 
chapter (2). 

4. Calculate the source terms for k and e 

5. Obtain the values of k and e at (n + 1) time level by solving equation 
(3.16) and (3.17) at the points excluding the near wall points. Where n 
refers to the psuedo time stepepping for k and e equtaion. For the near 
wall points apply wall-function treatment. 

6. Calculate the turbulent viscosity . 


7. update the boundary values of k and €. 
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8. Check for steady state by computing r.m.s value given by 

4*rms ~ 

if max(<j> rms ) < 10~ 3 goto step (2) else goto step (3). 



(3.45) 


In step (5) its important that the time step value (At) for k and e equation 
has to be smaller than that for the velocity so as to minimize the occurence 
of negative values of k and e (if some negative values of k and e very small in 
magnitude occurs then their absolute value is considered) . 


Chapter 4 


Results and Discussion 


The code is tested on two flow configurations, namely, two dimensional 
channel flow and the flow over a backward facing step. The results are com- 
pared with the published computational results. 

1. Channel Flow 

The flow in a channel has been simulated for two different Reynolds numbers, 
(i) Re=6600 and (n) Re=24000. The results due to Re=6600 are compared 
with the DNS results of Kim et al. 27 and those of the other Low-Reynolds 
number models of Chien 28 , Shih and Monsour 29 and Jones & Launder 30 . The 
results for Reynolds number of 24000 have been compared with the experi- 
mental data of Laufer 31 and the numerical results of Goldstein 33 and Lam k. 
Bremhorst 32 . 

For both the cases the Reynolds number is based on channel width. The 
flow configuration is shown in Figure (4.1). The boundary conditions of interest 
are 


at the inlet: 

U in = 1, V = 0 
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Figure 4.1: Turbulent flow through a channel 
k m = 1.5(L/j n I) 2 , where / = 0.05 

~3/4.3/2 

e in = ** A " , where A = 0.09 

9t = 0 

dx u 

at the outlet the gradient of the dependent variables in the flow direction is 
set to zero except for pressure, which has the conditionPout = 0. At the top 
and the bottom boundaries U = V = 0, and k and e the wall functions have 
been used 

For the case of Reynolds number of 6600, the distribution of turbulent 
kinetic energy (Fig. 4.2) and the nondimensional turbulent viscosity (Fig. 4.3) 
is compared with the results of DNS 27 and the other two low Reynolds number 
models 32 ’ 33 . Figure 4.2 depicts that the near-wall variation of turbulent kinetic 
energy from present computation corroborates well with the DNS data. How- 
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ever, the variation of turbulent viscosity predicted by the present computation 
does not match the DNS data so well (Fig. 4.3). For the case of Reynolds 
number of 24000, the near wall distribution of tangential velocity ( U + ) and ki- 
netic energy k from the prsent computation is compared with the experimental 
results of Laufer 31 and the numerical results of Cho and Goldstein 33 and Lam 
& Bremhorst 32 in Figures 4.4 and 4.5. The predicted results agree well with the 
experimental results of Laufer 31 and the other numerical results of Goldstein 33 . 

2 Backward-Facing Step flow 

Among the few selected flow situation that have been studied in detail is that 
of the flow over a backward facing step. It is often used as a test-case for flow 
solvers and turbulence models, since it embodies several difficult aspects of 
turbulent flows namely, flow seperation, reattachment, and the presence of the 
secondary recirculation regions. 

A schematic of the basic flow field is shown in Figure (4.6). The flow 
seperates at the step corner and reattaches at a distance Xr from the step. This 
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Figure 4.5: Near wall variation of turbulent kinetic energy for Re=24000 

distance is a function of the Reynolds number and the expansion ratio (^ i ). 
Several recent works (reference: 35-37) have studied the flow past a backward 
facing step and provide an excellent review of the current state of the art from 
both experimental and computational points of view. Among these, Speziale 34 
have applied a non-linear model for predicting the fluid flow and has found 
that introducing of anisotropy in the model leads to reduction in the error of 
predicting the reattachment length. The work presented here shows some of 
the features of standard k — e model as far the secondary recirculation zone is 
concerned. 

The fully developed flow past a backward facing step with an expansion 
ratio of 3:2, at a Reynold number of 46000 based on step height is chosen as a 
the primary reference. The computation are performed with a 110 x 73 finite 
volumes. The domain here is resolved finely because many earlier computations 
had serious error of around 25-30% due to insufficient resolution (Thangam 
and Hur 38 ). Regarding the boundary condition, a uniform velocity profile 
is prescribed at the inlet which is 11 H length upstream of the step. The 
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Figure 4.6: Flow past a backward facing step 

turbulent intensity I is 0.1, for e the inlet condition is specified as described 
in the channel flow problem. The downstream length is 24H so as to ensure 
fully developed condition. At the outlet boundary condition is as described 
in chapter3. No-slip boundary conditions for velocity and the wall function 
treatment for k and c are used at the walls. The results obtained for the velocity 
are shown in Figure (4.7) with a comparison with the experimental results of 
Kim et al . 36 , and the RNG results of Thangam and Speziale 37 . It was observed 
that a steady-state condition was not reached in the present computaions and 
that the flow remained guass-periodic as has been experimentally observed. 
Therefore the results shown in the figures are averages over a time cyle. 

The comparison of Reynold stresses (r IV ) is also made as shown in Figure 
(4.8) with the experimental results of Kim et al. 36 , second-order closure model 
calculation of Cellinligil and Mellor 35 , and the RNG calculation of Thangam 
and Speziale 38 . The Experimentas have not provided reliable data for Reynold 
stresses in the recirculation region and so a comparison in that region is not 



u 

o Experimental (Kim et ai. 36 ) 

. . . RNG (Speziale and Thangam 37 ) 
present 


Figure 4.7: Comparison of flow field 
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made. The RNG model predicted a reattachment length of 6.0£f while a second 
order model predicted a reattachment length of 7.7 H. and Thangam and Hur 
predicted the value of 5.58 H. The reattachment length found out by experi- 
ments is 7.0 H and the reattachment length obtained by present computation 
is 6.0H. 

The present computation seems to have been able in capturing the sec- 
ondary recirculation zone which other computations using the standard k — e 
model failed(references:38,39). The length of secondary recirculation is also 
predicted well as 0.83 H, as compared to 1.0 H obatined by Lien and Leschziner 39 
with a Reynold stress transport model. The variation of velocity (Figure 
4.9) near the bottom wall along the X direction shows the length of secondary 
and the primary recirculation zone. 

Also a comparison of turbulent viscosity with the experimental data of 
Driver and Segmiller 40 and the computational results done by the same au- 
thors using a two equation model of Jones and Launder is shown in Figure 
(4.10). The vector plots are shown in Figure (4.11-4.12). Figure (4.12)shows 
the secondary recirculation region. Figure (4.13 and 4.14) shows the stream- 
line plot for backwrad-facing step flow. Figure (4.14) also clearly shows the 
occurence of secondary recirculation. In fact a more fine grid may clearly show 
a tertiary eddy that seems to be nesting at the corner of the step. The com- 
putaional time taken by both the cases was 12.45 hours for channel flow and 
67.20 hours for backward-facing step flow on a DEC Alpha-3000 (version 2.2 ) 
dual processor machine. 

Conclusion 

A standard fc - e model is incorporated using a one layer near wall treatment 
in a Finite Volume code for flows in irregular geometries and is applied to two 
flow cases for validation with the following conclusions. 

1) The code is been able to predict most of the quantities in a channel flow in 


42 


a good agreement with the experimental and DNS results. 

2) The code is been able to solve the flow over a backward facing step in a 
good agreement with other published results. 

3) The computations seem to capture secondary recirculation regions better 
than other attempts with the standard k- e model. 

4) More investigation has to be done to find out the appropriate reason for 
secondary recirculation near the step corner. 
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> Experimental (Kim et a/. 36 ) 

< RNG (Speziale and Thangam 37 ) 

o Second-order closure model (Cellinligil and Mellor 35 ) 
present 


Figure 4.8: Comparison of turbulent shear stresses 








n 

v 


< Low Reynolds number k — e model (Jones and Launder 3 ) 
o Experiment (Driver and Seegmiller 40 ) 

present 


Figure 4.10: Distribution of eddy viscosity 
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